In [1]:
%matplotlib inline
%load_ext autoreload
%autoreload 2
import logging
logging.basicConfig(level=logging.ERROR)
logging.root.setLevel('INFO')

# Os
import pathlib
import sys
import gc
import os
import io
import json
import yaml
from marshmallow import pprint

# Data display
import qgrid
import IPython

# Scientific libs
import numpy as np
import pandas as pd
import scipy
import scipy.stats

# Geoprocissing libs
import geopandas as gpd
import rasterio as rio
import rasterio.crs
import fiona as fio
import fiona.crs
import shapely.geometry as sg

# Plotting
import matplotlib.pyplot as plt
from cycler import cycler

# EWS Library

assert os.path.isdir(os.environ["EWS_LIB_PATH"])
if not os.environ["EWS_LIB_PATH"] in sys.path:
    sys.path.append(os.environ["EWS_LIB_PATH"])    
from core.archive.git import get_last_commit
print(get_last_commit())
import core.utils.scitools as EWSSciTools
from core.docker import from_win_path_to_container as norm_pathname
from core.io.yaml import load_file, yaml, Dumper
from collections import namedtuple, OrderedDict, defaultdict
from core import init_notebook
from core.utils.ipython import display_end_of_script, nbconvert_html

paths = [os.path.join(os.environ['EWS_LIB_PATH'], x) for x in ["core", "wind"]]
import core.utils.io_utils
libSize = core.utils.io_utils.naturalsize(np.array([core.utils.io_utils.get_size(start_path=x) for x in paths]).sum())
print("EWS Library size: ", libSize)
del libSize, paths





def side_by_side(*args):
    #https://pandas.pydata.org/pandas-docs/stable/user_guide/style.html
    html_str = ''
    def process(x):
        if hasattr(x, "style"):
            return x.style.set_table_attributes("style='display:inline'")
        elif hasattr(x, "set_table_attributes"):
            return x.set_table_attributes("style='display:inline'")
        else:
            return x
    import pandas.io.formats.style
    for df in args:
        if isinstance(df, (pd.DataFrame, pd.io.formats.style.Styler)):
            html_str += process(df)._repr_html_()
        elif isinstance(df, tuple):
            html_str += process(df[0]).set_caption(df[1])._repr_html_()
    from core.utils.html import minify
    html_str = minify(html_str, engine="html_minify")
    IPython.display.display_html(html_str,raw=True)

def create_plot(**kws):
    from core.plot.eval_plot.eval_plot import EvalPlot
    from core.export.plot_base import PlotExporterBase
    logging.getLogger("core.plot.eval_plot.eval_plot").setLevel(logging.INFO)
    EWSTheme = PlotExporterBase.default_themes()["EWSTheme"]
    plot = EvalPlot(
        config={
            'pre_plot_hooks': ['{{ figure }}.add_subplot(111)'],
            #'despine': dict(enable=True, top=False, right=False, left=False, bottom=False),
            **kws
        },
        theme=EWSTheme).render()
    ax = plot.figure.axes[0]
    return plot, ax

init_notebook()
IPython.display.HTML('''<script>
code_show = true; 
function code_toggle() {
    if (code_show){
         $('div.input').hide();
    } else {
         $('div.input').show();
    }
    code_show = !code_show
}
</script>
<button type="button" onclick="code_toggle()"><b>TOGGLE CODE CELLS</b></button>''')
INFO:macropy.core.import_hooks:Expand macros in /usr/local/lib/python3.6/site-packages/macropy/core/hquotes.py
INFO:macropy.core.macros:Finding macros in 'macropy.core.hquotes'
INFO:macropy.core.macros:Importing macros from 'macropy.core.quotes' into 'macropy.core.hquotes'
INFO:macropy.core.import_hooks:Expand macros in /usr/local/lib/python3.6/site-packages/macropy/core/quotes.py
INFO:macropy.core.macros:Finding macros in 'macropy.core.quotes'
INFO:macropy.core.import_hooks:Expand macros in /usr/local/lib/python3.6/site-packages/macropy/core/failure.py
INFO:macropy.core.macros:Finding macros in 'macropy.core.failure'
INFO:macropy.core.macros:Importing macros from 'macropy.core.hquotes' into 'macropy.core.failure'
WARNING:core.eplot.patch:Patching Echarts for minification
WARNING:core.eplot.patch:Patching pandas for Eplot
WARNING:core.dataframe.patch_df.series:Patching pandas series for MultiIndex selection
WARNING:core.dataframe.patch_df.dataframe:Patching pandas dataframe for MultiIndex selection
WARNING:core.dataframe.patch_df.helpers:Patching pandas DatetimeIndex and Timestamp
WARNING:core.leaflet.patch:Patching geopandas for Leaflet rendering
commit 0199af867c24acea693109d8f7840b92f23b101d
Author: Farella Fabien <f.farella@ews-consulting.at>
Date:   Thu Nov 14 09:32:31 2019 +0100

    Added utility function to treat 3d lines/polygns


EWS Library size:  7.471 MB
Out[1]:
In [2]:
from core.docker import from_container_path_to_win
print(from_container_path_to_win(os.getcwd()))
T:\EWS_Projects\Gral_Gramm\askervein

BUILD PATHS

In [3]:
from core.utils.io_utils import make_path
WINDPRO_PATH = make_path(norm_pathname(r'T:\EWS_Projects\Gral_Gramm\askervein\windpro'))
if not WINDPRO_PATH.is_dir():
    WINDPRO_PATH.mkdir(parents=True, exist_ok=True)
BASIS_PATH = WINDPRO_PATH.joinpath("basis")
if not BASIS_PATH.is_dir():
    BASIS_PATH.mkdir(parents=True, exist_ok=True)

MAPS_PATH = BASIS_PATH.joinpath("karten")
if not MAPS_PATH.is_dir():
    MAPS_PATH.mkdir(parents=True, exist_ok=True)
DTM_PATH = BASIS_PATH.joinpath("höhenmodell")
if not DTM_PATH.is_dir():
    DTM_PATH.mkdir(parents=True, exist_ok=True)
ROUGH_PATH = BASIS_PATH.joinpath("rauigkeit")
if not ROUGH_PATH.is_dir():
    ROUGH_PATH.mkdir(parents=True, exist_ok=True)
COORDS_PATH = BASIS_PATH.joinpath("koordinaten")
if not COORDS_PATH.is_dir():
    COORDS_PATH.mkdir(parents=True, exist_ok=True)
list(BASIS_PATH.glob("*"))
Out[3]:
[PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/höhenmodell'),
 PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/karten'),
 PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/koordinaten'),
 PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/rauigkeit')]
In [11]:
import pathlib
DATA_PATH = pathlib.Path(norm_pathname(r"T:\EWS_Projects\Gral_Gramm\askervein\input_data"))
DATA_PATH.is_dir()
Out[11]:
True
In [23]:
import shapely.geometry
csv_fname = DATA_PATH.joinpath("Askervein_devc_loc.csv")
gdf = pd.read_csv(csv_fname)
gdf  =gdf[["Name", "AGL"]].assign(geometry=gdf.iloc[:, -2:].apply(shapely.geometry.Point, axis=1))
gdf = gpd.GeoDataFrame(gdf, crs=fio.crs.from_epsg(27700))
gdf.to_file("Askervein_locs.shp")
turbines_df = gdf.copy()
#gdf = gdf.to_crs(epsg=4326)
gdf
Out[23]:
Name AGL geometry
0 RS tower 1 m 1 POINT (74300 820980)
1 RS tower 3 m 3 POINT (74300 820980)
2 RS tower 5 m 5 POINT (74300 820980)
3 RS tower 8 m 8 POINT (74300 820980)
4 RS tower 15 m 15 POINT (74300 820980)
5 RS tower 24 m 24 POINT (74300 820980)
6 RS tower 34 m 34 POINT (74300 820980)
7 RS tower 49 m 49 POINT (74300 820980)
8 Cemetary BM -10000 POINT (73380 821940)
9 Milestone -10000 POINT (75284 821714)
10 TV tower -10000 POINT (73590 821711)
11 Base station -10000 POINT (74846 823306)
12 Pacman lake -10000 POINT (76207 823587)
13 Askernish House -10000 POINT (73654 823845)
14 HT Hill Top 10 POINT (75383 823737)
15 HT 10 m mf 10 POINT (75387 823735)
16 HT 10 m t 10 POINT (75381 823745)
17 HT tower 1 m 1 POINT (75381 823753)
18 HT tower 3 m 3 POINT (75381 823753)
19 HT tower 5 m 5 POINT (75381 823753)
20 HT tower 8 m 8 POINT (75381 823753)
21 HT tower 15 m 15 POINT (75381 823753)
22 HT tower 24 m 24 POINT (75381 823753)
23 HT tower 34 m 34 POINT (75381 823753)
24 HT tower 49 m 49 POINT (75381 823753)
25 HT WM 10 POINT (75360 823763)
26 CP Centre Point 10 POINT (75678 823465)
27 CP UK mf 10 POINT (75686 823457)
28 CP BSE 40 10 POINT (75680 823464)
29 CP FRG mf 10 POINT (75676 823466)
... ... ... ...
44 BSE130 10 POINT (76338 822860)
45 BSE150 10 POINT (76490 822723)
46 BSE170 10 POINT (76636 822585)
47 ANE10 10 POINT (75454 823812)
48 ANE20 10 POINT (75523 823884)
49 ANE40 10 POINT (75661 824017)
50 ASW10 10 POINT (75319 823667)
51 ASW20 10 POINT (75257 823604)
52 ASW35 10 POINT (75162 823498)
53 ASW50 10 POINT (75050 823378)
54 ASW UK 30 m tower 30 POINT (74970 823295)
55 ASW85 10 POINT (74813 823122)
56 AANE10 10 POINT (75746 823540)
57 AANE20 10 POINT (75807 823610)
58 AANE30 10 POINT (75871 823675)
59 AANE40 10 POINT (75938 823745)
60 AANE50 10 POINT (76006 823813)
61 AANE60 10 POINT (76073 823886)
62 AASW10 mf 10 POINT (75613 823396)
63 AASW10 t 10 POINT (75623 823388)
64 AASW20 10 POINT (75553 823320)
65 AASW30 mf 10 POINT (75483 823253)
66 AASW30 t 10 POINT (75493 823244)
67 AASW40 10 POINT (75417 823174)
68 AASW50 mf 10 POINT (75341 823107)
69 AASW50 t 10 POINT (75352 823100)
70 AASW60 10 POINT (75274 823038)
71 AASW70 10 POINT (75208 822968)
72 AASW80 10 POINT (75140 822895)
73 AASW90 10 POINT (75069 822820)

74 rows × 3 columns

In [24]:
gdf.to_leaflet(display=True);
In [25]:
center_pt = gdf.geometry.buffer(25000).unary_union.centroid
list(map(int, (center_pt.x, center_pt.y)))
Out[25]:
[74865, 822670]

BUILD CONFIGURATION

In [54]:
EPSG_CODE = 27700
ANALYSIS_CFG_str = rf'''

projection:
    epsg: {EPSG_CODE}
    proj4s: {rio.crs.CRS.from_epsg(EPSG_CODE).to_proj4()}
    crs: {rio.crs.CRS.from_epsg(EPSG_CODE).to_dict()!r}
    name: {rio.crs.CRS.from_epsg(EPSG_CODE).to_wkt().split('"')[1].split('"')[0]}

domains:
    items:
        - domain_name: small_domain
          buffer_distance: 5000
          dtm_resolution: 5
          
        - domain_name: large_domain
          buffer_distance: 10000
          dtm_resolution: 50
          
dtm:
    path: {str(DTM_PATH)}
    items:
        - domain_name: large_domain
          parts:
              - name: SRTM
                filename:  S:\SRTM_Data\all_srtms.vrt
                
        - domain_name: small_domain
          parts:
              - name: SRTM
                filename:  S:\SRTM_Data\all_srtms.vrt
          
roughness:
    path: {str(ROUGH_PATH)}
    items:
        - domain_name: large_domain
        - domain_name: small_domain

maps:   
    path: {str(MAPS_PATH)}
    items:
        - domain_name: large_domain
          mpp: 15
          server_name: 'Google Satellite'
          filename: 'google_satellite_large'
          
        - domain_name: large_domain
          mpp: 15
          server_name: 'ARCGIS World_Shaded_Relief'
          filename: 'shaded_relief_large'
          
        - domain_name: small_domain
          mpp: 8
          server_name: 'Google Satellite'
          filename: 'google_satellite_small' 

        
'''


ANALYSIS_CFG = yaml.safe_load(io.StringIO(ANALYSIS_CFG_str))
ANALYSIS_CFG["projection"] = namedtuple("Configuration", ANALYSIS_CFG["projection"].keys())(*ANALYSIS_CFG["projection"].values())
ANALYSIS_CFG = namedtuple("Configuration", ANALYSIS_CFG.keys())(*ANALYSIS_CFG.values())


print(ANALYSIS_CFG_str)


gdf = turbines_df.to_crs(ANALYSIS_CFG.projection.crs).set_index("Name").copy()
geom = gdf.geometry
domains_dict = dict()
for domain_cfg in ANALYSIS_CFG.domains.get("items", []):
    domain_name = domain_cfg["domain_name"]
    print(domain_name)
    
    r = float(domain_cfg["dtm_resolution"])
    buffer_distance = float(domain_cfg["buffer_distance"])
    excluded_names = domain_cfg.get("excluded_names", [])
    
    geom = gdf.drop(excluded_names, errors="ignore", axis=0).geometry
    b = np.array(geom.buffer(buffer_distance + 50000).unary_union.buffer(-50000).bounds)    
    download_shape = sg.box(*np.column_stack((r * np.floor(b[:2] / r), r * np.ceil(b[2:] / r))).T.flatten())
    gdf.loc[domain_name, "Type"] = "Domain" 
    gdf.loc[domain_name, "geometry"] = download_shape
    
    domains_dict[domain_name] = {
        "buffer_distance": buffer_distance,
        "dtm_resolution": r,
        "shape": download_shape,
        "bounds": download_shape.bounds,
    }
del geom, b, r, domain_name, buffer_distance, download_shape

gdf.reset_index().to_leaflet(buffer=100)[0]

projection:
    epsg: 27700
    proj4s: +init=epsg:27700
    crs: {'init': 'epsg:27700'}
    name: OSGB 1936 / British National Grid

domains:
    items:
        - domain_name: small_domain
          buffer_distance: 5000
          dtm_resolution: 5
          
        - domain_name: large_domain
          buffer_distance: 10000
          dtm_resolution: 50
          
dtm:
    path: /home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/höhenmodell
    items:
        - domain_name: large_domain
          parts:
              - name: SRTM
                filename:  S:\SRTM_Data\all_srtms.vrt
                
        - domain_name: small_domain
          parts:
              - name: SRTM
                filename:  S:\SRTM_Data\all_srtms.vrt
          
roughness:
    path: /home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/rauigkeit
    items:
        - domain_name: large_domain
        - domain_name: small_domain

maps:   
    path: /home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/karten
    items:
        - domain_name: large_domain
          mpp: 15
          server_name: 'Google Satellite'
          filename: 'google_satellite_large'
          
        - domain_name: large_domain
          mpp: 15
          server_name: 'ARCGIS World_Shaded_Relief'
          filename: 'shaded_relief_large'
          
        - domain_name: small_domain
          mpp: 8
          server_name: 'Google Satellite'
          filename: 'google_satellite_small' 

        

small_domain
large_domain
Out[54]:

DOWNLOAD DTM

In [42]:
from core.gis.raster.raster_base import RasterBase
import tempfile
import core.gis.raster.ops as RasterOps
from core.utils.system import execute_in_terminal
import shutil
dtm_dict = dict()
for r_config in  ANALYSIS_CFG.dtm.get("items", []):
    domain_name = r_config["domain_name"]
    domain_cfg = domains_dict[domain_name]
    resolution = domain_cfg["dtm_resolution"]
    bounds = domain_cfg["bounds"]
    raster_config = {
        "name": domain_name,
        "resampling": {
            "active": True,
            "method": "bilinear",
            "resolution": resolution,
            "bounds": list(bounds),
            "projection": {
              "proj4s": ANALYSIS_CFG.projection.proj4s,
              "name": ANALYSIS_CFG.projection.name
            }
        },
        "parts": r_config["parts"]
    }
    pprint(raster_config)    
    raster = RasterBase(config=raster_config)
    raster.read()
    dtm_dict[domain_name] = raster
    plt.imshow(raster.data, extent=raster.extent) 
    
    path = make_path(ANALYSIS_CFG.dtm.get("path", os.getcwd()))
    dtm_fname = path.joinpath(f"{domain_name}.tif")
    grd_fname = path.joinpath(f"{domain_name}.asc.txt")
    with tempfile.NamedTemporaryFile(delete=False) as t_file:
        RasterOps.data_to_geotiff(
            t_file.name, raster.data.astype("float32"), raster.transform, 
            profile={"dtype": "float32", "crs": raster.crs, "nodata":-999}
        )        
        shutil.move(t_file.name, dtm_fname)
    assert dtm_fname.exists()
    execute_in_terminal(f"gdal_translate -of AAIGrid '{dtm_fname}' '{grd_fname}'", verbose=True)
    #dtm_fname.unlink()
INFO:core.gis.raster.raster_base:Reading with resampling
INFO:core.gis.raster.raster_base:[58400.0, 806000.0, 91600.0, 839000.0]
INFO:core.gis.raster.raster_base:{'driver': 'VRT', 'dtype': 'int16', 'nodata': -32768.0, 'width': 664, 'height': 660, 'count': 1, 'crs': CRS.from_dict(init='epsg:27700'), 'transform': Affine(50.0, 0.0, 58400.0,
       0.0, -50.0, 839000.0), 'blockxsize': 512, 'blockysize': 128, 'tiled': True}
{'name': 'large_domain',
 'parts': [{'filename': 'S:\\SRTM_Data\\all_srtms.vrt', 'name': 'SRTM'}],
 'resampling': {'active': True,
                'bounds': [58400.0, 806000.0, 91600.0, 839000.0],
                'method': 'bilinear',
                'projection': {'name': 'OSGB 1936 / British National Grid',
                               'proj4s': '+init=epsg:27700'},
                'resolution': 50.0}}
gdal_translate -of AAIGrid '/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/höhenmodell/large_domain.tif' '/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/höhenmodell/large_domain.asc.txt'
Child returned 0
INFO:core.gis.raster.raster_base:Reading with resampling
INFO:core.gis.raster.raster_base:[68440.0, 816040.0, 81580.0, 828960.0]
0
{'name': 'small_domain',
 'parts': [{'filename': 'S:\\SRTM_Data\\all_srtms.vrt', 'name': 'SRTM'}],
 'resampling': {'active': True,
                'bounds': [68440.0, 816040.0, 81580.0, 828960.0],
                'method': 'bilinear',
                'projection': {'name': 'OSGB 1936 / British National Grid',
                               'proj4s': '+init=epsg:27700'},
                'resolution': 5.0}}
INFO:core.gis.raster.raster_base:{'driver': 'VRT', 'dtype': 'int16', 'nodata': -32768.0, 'width': 2628, 'height': 2584, 'count': 1, 'crs': CRS.from_dict(init='epsg:27700'), 'transform': Affine(5.0, 0.0, 68440.0,
       0.0, -5.0, 828960.0), 'blockxsize': 512, 'blockysize': 128, 'tiled': True}
gdal_translate -of AAIGrid '/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/höhenmodell/small_domain.tif' '/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/höhenmodell/small_domain.asc.txt'
Child returned 0
0
In [43]:
fig, axL = plt.subplots(1, len(ANALYSIS_CFG.dtm.get("items", [])), figsize=(16, 12) )
axL = np.array(axL).flatten()
rdict = dict()
for i, r_config in  enumerate(ANALYSIS_CFG.dtm.get("items", [])):
    ax = axL[i]
    domain_name = r_config["domain_name"]
    path = make_path(ANALYSIS_CFG.dtm.get("path", os.getcwd()))
    dtm_fname = path.joinpath(f"{domain_name}.tif")
    with rio.open(str(dtm_fname)) as src:
        data = src.read(1)
        bounds = src.bounds
        extent = src.bounds[::2] + src.bounds[1::2]
    cmap = "gist_earth"
    norm = None
    rdict[domain_name] = data
    im = ax.imshow(data, cmap=cmap, norm=norm, extent=extent)
    plt.colorbar(im, ax=ax, orientation="horizontal")    
    if True:
        ax.plot(
            [pt.x for pt in gdf.query("Type != 'Domain'").geometry], 
            [pt.y for pt in gdf.query("Type != 'Domain'").geometry], 
            mfc="w", color="k", marker="o", ls="none"
        )
    ax.set(xlim=extent[:2], ylim=extent[2:], aspect="equal", title=domain_name)
fig.tight_layout();

DOWNLOAD CORINE LANDCOVER

In [55]:
import rasterio.features
import requests, urllib3
import shutil
from core.utils.system import execute_in_terminal
In [56]:
CLASS_LUT = r"""ID	Land Cover Type	EMD Roughness Length
1	Continuous urban fabric 1.1.1	z0=0.50000
2	Discontinuous urban fabric 1.1.2	z0=0.40000
3	Industrial or commercial units 1.2.1	z0=0.70000
4	Road and rail networks and associated land 1.2.2	z0=0.10000
5	Port areas 1.2.3	z0=0.50000
6	Airports 1.2.6	z0=0.03000
7	Mineral extraction sites 1.3.1	z0=0.10000
8	Dump sites 1.3.2	z0=0.10000
9	Construction sites 1.3.3	z0=0.30000
10	Green urban areas 1.4.1	z0=0.40000
11	Sport and leisure facilities 1.4.2	z0=0.50000
12	Non-irrigated arable land 2.1.1	z0=0.05600
13	Permanently irrigated land 2.1.2	z0=0.05600
14	Rice fields 2.1.3	z0=0.01840
15	Vineyards 2.2.1	z0=0.30000
16	Fruit trees and berry plantations 2.2.2	z0=0.40000
17	Olive groves 2.2.3	z0=0.40000
18	Pastures 2.3.1	z0=0.03600
19	Annual crops associated with permanent crops 2.4.1	z0=0.05600
20	Complex cultivation patterns 2.4.2	z0=0.05600
21	Land principally occupied by agriculture, with significant areas of natural vegetation 2.4.3	z0=0.05600
22	Agro-forestry areas 2.4.4	z0=0.50000
23	Broad-leaved forest 3.1.1	z0=0.50000
24	Coniferous forest 3.1.2	z0=0.50000
25	Mixed forest 3.1.3	z0=0.50000
26	Natural grassland 3.2.1	z0=0.05600
27	Moors and heathland 3.2.2	z0=0.06000
28	Sclerophyllous vegetation 3.2.3	z0=0.05600
29	Sparse coniferous forest, cc. 10-30% 3.2.4	z0=0.40000
30	Beaches, dunes, and sand plains 3.3.1	z0=0.01000
31	Bare rock 3.3.2	z0=0.05000
32	Sparsely vegetated areas 3.3.3	z0=0.20000
33	Burnt areas 3.3.4	z0=0.20000
34	Glaciers and perpetual snow 3.3.5	z0=0.20000
35	Inland marshes 4.1.1	z0=0.05000
36	Peat bogs 4.1.2	z0=0.01840
37	Salt marshes 4.2.1	z0=0.03480
38	Salines 4.2.2	z0=0.03000
39	Intertidal flats 4.2.3	z0=0.00050
40	River 5.1.1	z0=0.00000
41	Lake 5.1.2	z0=0.00000
42	Coastal lagoons 5.2.1	z0=0.00000
43	Estuaries 5.2.2	z0=0.00000
44	Sea and ocean 5.2.3	z0=0.00000
45	NO DATA	z0=-1.00000"""
CLASS_LUT = pd.read_csv(io.StringIO(CLASS_LUT), sep="\t")
CLASS_LUT['EMD Roughness Length'] = [_[1] for _ in CLASS_LUT['EMD Roughness Length'].str.split("z0=")]
CLASS_LUT['z0'] = CLASS_LUT['EMD Roughness Length'].astype("float")
del CLASS_LUT['EMD Roughness Length']
CLASS_LUT["ID"] = CLASS_LUT["ID"].astype("int")

CLASS_LUT["Code_18"] = [_[-1] for _ in CLASS_LUT['Land Cover Type'].str.split(" ")]
CLASS_LUT["Code_18"] = CLASS_LUT["Code_18"].str.replace(".", "").replace("DATA", "-1").astype("int")


CLASS_LUT.idisplay;
CLASS_LUT_DF = CLASS_LUT.set_index("Code_18")
CLASS_LUT = CLASS_LUT.set_index("Code_18")["z0"].to_dict()
label_dict = CLASS_LUT_DF.apply(
    lambda row: f"{row['Land Cover Type']} - z0={row.z0}", axis=1
).to_dict()
label_dict
ID Land Cover Type z0 Code_18
0 1 Continuous urban fabric 1.1.1 0.5000 111
1 2 Discontinuous urban fabric 1.1.2 0.4000 112
2 3 Industrial or commercial units 1.2.1 0.7000 121
3 4 Road and rail networks and associated land 1.2.2 0.1000 122
4 5 Port areas 1.2.3 0.5000 123
5 6 Airports 1.2.6 0.0300 126
6 7 Mineral extraction sites 1.3.1 0.1000 131
7 8 Dump sites 1.3.2 0.1000 132
8 9 Construction sites 1.3.3 0.3000 133
9 10 Green urban areas 1.4.1 0.4000 141
10 11 Sport and leisure facilities 1.4.2 0.5000 142
11 12 Non-irrigated arable land 2.1.1 0.0560 211
12 13 Permanently irrigated land 2.1.2 0.0560 212
13 14 Rice fields 2.1.3 0.0184 213
14 15 Vineyards 2.2.1 0.3000 221
15 16 Fruit trees and berry plantations 2.2.2 0.4000 222
16 17 Olive groves 2.2.3 0.4000 223
17 18 Pastures 2.3.1 0.0360 231
18 19 Annual crops associated with permanent crops 2... 0.0560 241
19 20 Complex cultivation patterns 2.4.2 0.0560 242
20 21 Land principally occupied by agriculture, with... 0.0560 243
21 22 Agro-forestry areas 2.4.4 0.5000 244
22 23 Broad-leaved forest 3.1.1 0.5000 311
23 24 Coniferous forest 3.1.2 0.5000 312
24 25 Mixed forest 3.1.3 0.5000 313
25 26 Natural grassland 3.2.1 0.0560 321
26 27 Moors and heathland 3.2.2 0.0600 322
27 28 Sclerophyllous vegetation 3.2.3 0.0560 323
28 29 Sparse coniferous forest, cc. 10-30% 3.2.4 0.4000 324
29 30 Beaches, dunes, and sand plains 3.3.1 0.0100 331
30 31 Bare rock 3.3.2 0.0500 332
31 32 Sparsely vegetated areas 3.3.3 0.2000 333
32 33 Burnt areas 3.3.4 0.2000 334
33 34 Glaciers and perpetual snow 3.3.5 0.2000 335
34 35 Inland marshes 4.1.1 0.0500 411
35 36 Peat bogs 4.1.2 0.0184 412
36 37 Salt marshes 4.2.1 0.0348 421
37 38 Salines 4.2.2 0.0300 422
38 39 Intertidal flats 4.2.3 0.0005 423
39 40 River 5.1.1 0.0000 511
40 41 Lake 5.1.2 0.0000 512
41 42 Coastal lagoons 5.2.1 0.0000 521
42 43 Estuaries 5.2.2 0.0000 522
43 44 Sea and ocean 5.2.3 0.0000 523
44 45 NO DATA -1.0000 -1
Out[56]:
{111: 'Continuous urban fabric 1.1.1 - z0=0.5',
 112: 'Discontinuous urban fabric 1.1.2 - z0=0.4',
 121: 'Industrial or commercial units 1.2.1 - z0=0.7',
 122: 'Road and rail networks and associated land 1.2.2 - z0=0.1',
 123: 'Port areas 1.2.3 - z0=0.5',
 126: 'Airports 1.2.6 - z0=0.03',
 131: 'Mineral extraction sites 1.3.1 - z0=0.1',
 132: 'Dump sites 1.3.2 - z0=0.1',
 133: 'Construction sites 1.3.3 - z0=0.3',
 141: 'Green urban areas 1.4.1 - z0=0.4',
 142: 'Sport and leisure facilities 1.4.2 - z0=0.5',
 211: 'Non-irrigated arable land 2.1.1 - z0=0.056',
 212: 'Permanently irrigated land 2.1.2 - z0=0.056',
 213: 'Rice fields 2.1.3 - z0=0.0184',
 221: 'Vineyards 2.2.1 - z0=0.3',
 222: 'Fruit trees and berry plantations 2.2.2 - z0=0.4',
 223: 'Olive groves 2.2.3 - z0=0.4',
 231: 'Pastures 2.3.1 - z0=0.036',
 241: 'Annual crops associated with permanent crops 2.4.1 - z0=0.056',
 242: 'Complex cultivation patterns 2.4.2 - z0=0.056',
 243: 'Land principally occupied by agriculture, with significant areas of natural vegetation 2.4.3 - z0=0.056',
 244: 'Agro-forestry areas 2.4.4 - z0=0.5',
 311: 'Broad-leaved forest 3.1.1 - z0=0.5',
 312: 'Coniferous forest 3.1.2 - z0=0.5',
 313: 'Mixed forest 3.1.3 - z0=0.5',
 321: 'Natural grassland 3.2.1 - z0=0.056',
 322: 'Moors and heathland 3.2.2 - z0=0.06',
 323: 'Sclerophyllous vegetation 3.2.3 - z0=0.056',
 324: 'Sparse coniferous forest, cc. 10-30% 3.2.4 - z0=0.4',
 331: 'Beaches, dunes, and sand plains 3.3.1 - z0=0.01',
 332: 'Bare rock 3.3.2 - z0=0.05',
 333: 'Sparsely vegetated areas 3.3.3 - z0=0.2',
 334: 'Burnt areas 3.3.4 - z0=0.2',
 335: 'Glaciers and perpetual snow 3.3.5 - z0=0.2',
 411: 'Inland marshes 4.1.1 - z0=0.05',
 412: 'Peat bogs 4.1.2 - z0=0.0184',
 421: 'Salt marshes 4.2.1 - z0=0.0348',
 422: 'Salines 4.2.2 - z0=0.03',
 423: 'Intertidal flats 4.2.3 - z0=0.0005',
 511: 'River 5.1.1 - z0=0.0',
 512: 'Lake 5.1.2 - z0=0.0',
 521: 'Coastal lagoons 5.2.1 - z0=0.0',
 522: 'Estuaries 5.2.2 - z0=0.0',
 523: 'Sea and ocean 5.2.3 - z0=0.0',
 -1: 'NO DATA - z0=-1.0'}
In [57]:
for r_config in  ANALYSIS_CFG.roughness.get("items", []):
    domain_name = r_config["domain_name"]
    domain_cfg = domains_dict[domain_name]
    #resolution = domain_cfg["dtm_resolution"]
    #bounds = domain_cfg["bounds"]
    raster = dtm_dict[domain_name]
    resolution = raster.resolution[0]
    bounds = raster.bounds
    raster.read()

    
    path = make_path(ANALYSIS_CFG.roughness.get("path", os.getcwd()))
    
    base_url = "https://image.discomap.eea.europa.eu/arcgis/rest/services/Corine/CLC2018_WM/MapServer/0/query?"

    base_kws = {
        "geometry": ",".join(map(str, bounds)),
        "geometryType": "esriGeometryEnvelope",
        "inSR": ANALYSIS_CFG.projection.epsg,
        "spatialRel": "esriSpatialRelIntersects",    
        "returnTrueCurves": True,
        "outSR":  ANALYSIS_CFG.projection.epsg,
        "returnZ": False,
        "returnM": False,
        "f": "json",
        "text":  "",
        "time":  "",
        "relationParam":  "",
        "outFields":  "",
        "maxAllowableOffset":  "",
        "geometryPrecision":  "",
        "orderByFields":  "",
        "groupByFieldsForStatistics":  "",
        "outStatistics":  "",
        "gdbVersion": "",
        "returnDistinctValues": False,
        "resultOffset":  "",
        "resultRecordCount":  "",
        "queryByDistance":  "",
        "returnExtentsOnly": False,
        "datumTransformation":  "",
        "parameterValues":  "",
        "rangeValues":  "",
    }
    url = base_url +  urllib3.request.urlencode({
        "returnIdsOnly": False,
        "returnCountOnly": True,
        "returnGeometry": False,
        **base_kws
    })
    #print(url)

    n_expected = requests.api.get(url).json()["count"]
    print(n_expected)
    url = base_url +  urllib3.request.urlencode({
        "returnIdsOnly": True,
        "returnCountOnly": False,
        "returnGeometry": False,
        **base_kws
    })

    ids = requests.api.get(url).json()["objectIds"]
    print(len(ids))
    assert len(ids) == n_expected, "Number of ids doesn't match!!!"
    MAX_RECORDS_COUNT = 1000
    n, r = divmod(len(ids), MAX_RECORDS_COUNT)
    n = n+1 if r > 0 else n
    dl_list = []
    for cpt in range(n):
        dl_list += [ids[cpt*MAX_RECORDS_COUNT:(cpt+1)*MAX_RECORDS_COUNT]]
    print([len(x) for x in dl_list])

    c_gdf = list()
    for idList in dl_list:
        s = ",".join(map(str, idList))
        url = base_url +  urllib3.request.urlencode({
            #"where": f"OBJECTID in ({s})",
            "where": f"(OBJECTID >= {idList[0]}) AND (OBJECTID <= {idList[-1]})",
            #"objectIds": ",".join(map(str, idList)),
            "returnIdsOnly": False,
            "returnCountOnly": False,
            "returnGeometry": True,
            **base_kws
        })
        print(url)
        c_gdf += [gpd.read_file(url)]
    len(c_gdf)
    c_gdf = gpd.GeoDataFrame(pd.concat(c_gdf, axis=0).astype({"Code_18": "int"}), crs= ANALYSIS_CFG.projection.crs) 
    c_gdf.info()
    

    assert c_gdf.index.size == n_expected, "Number of ids doesn't match!!!"
    
    c_gdf["geometry"] = c_gdf.geometry.intersection(domain_cfg["shape"].buffer(4*resolution))
    c_gdf = c_gdf[c_gdf.geometry.area > 0]
    
    
    c_gdf["RoughnessL"] = c_gdf.Code_18.astype("int").map(CLASS_LUT.get).replace([None], np.nan)
    c_gdf["Label"] = c_gdf.Code_18.astype("int").map(label_dict.get).replace([None], np.nan)
    print(c_gdf["RoughnessL"].isnull().sum())
    c_gdf = c_gdf[c_gdf.RoughnessL >= 0][["RoughnessL", "geometry", "Code_18", "Label"]]

    import requests
    legend_json = requests.api.get("https://image.discomap.eea.europa.eu/arcgis/rest/services/Corine/CLC2018_WM/MapServer/0?f=pjson").json()
    legend = pd.DataFrame.from_records(legend_json['drawingInfo']['renderer']['uniqueValueInfos'])[["value", "label"]].astype({"value": "int"}).set_index("value")["label"].to_dict()
    cdict = {int(v["value"]): tuple(np.array(v['symbol']["color"]) / 255.0)  for v in legend_json['drawingInfo']['renderer']['uniqueValueInfos']}
    legend
    def style_function(feature):
        import matplotlib as mpl
        c = mpl.colors.rgb2hex(cdict[feature['properties']['Code_18']][:3])
        return {'color': c, 'fillColor': c, 'fillOpacity': 0.4}


    fig, ax = plt.subplots(1, 1, figsize=(12, 12) )
    ax = c_gdf.plot(ax=ax, column="Code_18")
    ax.set(xlim=bounds[::2], ylim=bounds[1::2], aspect="equal")
    c_raster = rio.features.rasterize(
        c_gdf[["geometry", "Code_18"]].apply(tuple, axis=1),
        out_shape=raster.data.shape,
        transform=raster.transform,
        fill=1,
        all_touched=True,
        dtype=np.int32
    )
    fig, ax = plt.subplots(1, 1, figsize=(12, 12) )
    im = ax.imshow(c_raster, extent=raster.extent)
    ax.set(xlim=bounds[::2], ylim=bounds[1::2], aspect="equal")

    shp_fname = path.joinpath(f"{domain_name}_landcover_clc2018.shp")
    c_gdf.to_file(shp_fname, driver="ESRI Shapefile")
    clc_fname = path.joinpath(f"{domain_name}_landcover_clc2018.tif")
    grd_fname = path.joinpath(f"{domain_name}_landcover_clc2018.asc.txt")
    
    m = c_gdf.to_leaflet(style_function=style_function)[0]
    m.save(str(path.joinpath(f"{domain_name}_landcover_clc2018.html")))
    

    with tempfile.NamedTemporaryFile(delete=False) as t_file:
            RasterOps.data_to_geotiff(
                t_file.name, np.ma.masked_less(c_raster, 0).astype("int32"), raster.transform, 
                profile={"dtype": "int32", "crs": raster.crs, "nodata":-999}
            )            
            shutil.move(t_file.name, clc_fname)
    assert clc_fname.exists()    
    execute_in_terminal(f"gdal_translate -of AAIGrid '{clc_fname}' '{grd_fname}'", verbose=True)
INFO:core.gis.raster.raster_base:Reading with resampling
INFO:core.gis.raster.raster_base:[58400.0, 806000.0, 91600.0, 839000.0]
INFO:core.gis.raster.raster_base:{'driver': 'VRT', 'dtype': 'int16', 'nodata': -32768.0, 'width': 664, 'height': 660, 'count': 1, 'crs': CRS.from_dict(init='epsg:27700'), 'transform': Affine(50.0, 0.0, 58400.0,
       0.0, -50.0, 839000.0), 'blockxsize': 512, 'blockysize': 128, 'tiled': True}
91
91
[91]
https://image.discomap.eea.europa.eu/arcgis/rest/services/Corine/CLC2018_WM/MapServer/0/query?where=%28OBJECTID+%3E%3D+1626174%29+AND+%28OBJECTID+%3C%3D+2372333%29&returnIdsOnly=False&returnCountOnly=False&returnGeometry=True&geometry=58400.0%2C806000.0%2C91600.0%2C839000.0&geometryType=esriGeometryEnvelope&inSR=27700&spatialRel=esriSpatialRelIntersects&returnTrueCurves=True&outSR=27700&returnZ=False&returnM=False&f=json&text=&time=&relationParam=&outFields=&maxAllowableOffset=&geometryPrecision=&orderByFields=&groupByFieldsForStatistics=&outStatistics=&gdbVersion=&returnDistinctValues=False&resultOffset=&resultRecordCount=&queryByDistance=&returnExtentsOnly=False&datumTransformation=&parameterValues=&rangeValues=
<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 91 entries, 0 to 90
Data columns (total 2 columns):
Code_18     91 non-null int64
geometry    91 non-null object
dtypes: int64(1), object(1)
memory usage: 1.5+ KB
0
gdal_translate -of AAIGrid '/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/rauigkeit/large_domain_landcover_clc2018.tif' '/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/rauigkeit/large_domain_landcover_clc2018.asc.txt'
Child returned 0
INFO:core.gis.raster.raster_base:Reading with resampling
INFO:core.gis.raster.raster_base:[68440.0, 816040.0, 81580.0, 828960.0]
INFO:core.gis.raster.raster_base:{'driver': 'VRT', 'dtype': 'int16', 'nodata': -32768.0, 'width': 2628, 'height': 2584, 'count': 1, 'crs': CRS.from_dict(init='epsg:27700'), 'transform': Affine(5.0, 0.0, 68440.0,
       0.0, -5.0, 828960.0), 'blockxsize': 512, 'blockysize': 128, 'tiled': True}
0
41
41
[41]
https://image.discomap.eea.europa.eu/arcgis/rest/services/Corine/CLC2018_WM/MapServer/0/query?where=%28OBJECTID+%3E%3D+1627310%29+AND+%28OBJECTID+%3C%3D+2372333%29&returnIdsOnly=False&returnCountOnly=False&returnGeometry=True&geometry=68440.0%2C816040.0%2C81580.0%2C828960.0&geometryType=esriGeometryEnvelope&inSR=27700&spatialRel=esriSpatialRelIntersects&returnTrueCurves=True&outSR=27700&returnZ=False&returnM=False&f=json&text=&time=&relationParam=&outFields=&maxAllowableOffset=&geometryPrecision=&orderByFields=&groupByFieldsForStatistics=&outStatistics=&gdbVersion=&returnDistinctValues=False&resultOffset=&resultRecordCount=&queryByDistance=&returnExtentsOnly=False&datumTransformation=&parameterValues=&rangeValues=
<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 41 entries, 0 to 40
Data columns (total 2 columns):
Code_18     41 non-null int64
geometry    41 non-null object
dtypes: int64(1), object(1)
memory usage: 736.0+ bytes
0
gdal_translate -of AAIGrid '/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/rauigkeit/small_domain_landcover_clc2018.tif' '/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/rauigkeit/small_domain_landcover_clc2018.asc.txt'
Child returned 0
0
In [59]:
def get_corine_cmap():
    import matplotlib as mpl
    res = list()
    for cfg in requests.get(
        "https://image.discomap.eea.europa.eu/arcgis/rest/services/Corine/CLC2018_WM/MapServer/0?f=pjson"
    ).json()["drawingInfo"]["renderer"]["uniqueValueInfos"]:
        value = cfg["value"]
        label = cfg["label"]
        color = cfg["symbol"]["color"][:-1]
        res.append((label, value, color))
    res = pd.DataFrame(
        res, columns=["Label", "Code", "RGBColor"]).set_index("Code")
    res["HEXColor"] = res.RGBColor.apply(
        lambda ls: mpl.colors.to_hex(np.array(ls) / 255)).str.upper()
    
    cmap = mpl.colors.ListedColormap(res["HEXColor"].tolist())
    norm = mpl.colors.BoundaryNorm(pd.Float64Index(res.index).tolist(), cmap.N) 
        
    return cmap, norm


corine_cmap, corine_norm = get_corine_cmap()
fig, axL = plt.subplots(1, len(ANALYSIS_CFG.roughness.get("items", [])), figsize=(16, 12) )
axL = np.array(axL).flatten()
for i, r_config in  enumerate(ANALYSIS_CFG.roughness.get("items", [])):
    ax = axL[i]
    domain_name = r_config["domain_name"]
    path = make_path(ANALYSIS_CFG.roughness.get("path", os.getcwd()))
    dtm_fname = path.joinpath(f"{domain_name}_landcover_clc2018.asc.txt")
    with rio.open(str(dtm_fname)) as src:
        data = src.read(1)
        bounds = src.bounds
        extent = src.bounds[::2] + src.bounds[1::2]
    im = ax.imshow(data, cmap=corine_cmap, norm=corine_norm, extent=extent)
    plt.colorbar(im, ax=ax, orientation="horizontal")    
    if True:
        ax.plot(
            [pt.x for pt in gdf.query("Type != 'Domain'").geometry], 
            [pt.y for pt in gdf.query("Type != 'Domain'").geometry], 
            mfc="w", color="k", marker="o", ls="none"
        )
    ax.set(xlim=extent[:2], ylim=extent[2:], aspect="equal", title=domain_name)
fig.tight_layout();
In [60]:
import core.gis.raster.ops as RasterOps
idx = gdf.query("Type != 'Domain'").index
gdf.loc[idx, "TerrainHeight"] = RasterOps.interpolate_on_data(
    gdf.loc[idx].geometry.apply(lambda pt: pd.Series([pt.x, pt.y])),
    dtm_dict["small_domain"].data,
    dtm_dict["small_domain"].transform,
    method="linear",
    bounds_error=False,
    fill_value=np.nan
).round(1)
gdf.TerrainHeight
Out[60]:
Name
RS tower 1 m           0.0
RS tower 3 m           0.0
RS tower 5 m           0.0
RS tower 8 m           0.0
RS tower 15 m          0.0
RS tower 24 m          0.0
RS tower 34 m          0.0
RS tower 49 m          0.0
Cemetary BM           14.2
Milestone              7.0
TV tower               6.3
Base station           9.0
Pacman lake           38.0
Askernish House        4.0
HT Hill Top          122.0
HT 10 m mf           122.0
HT 10 m t            122.0
HT tower 1 m         122.0
HT tower 3 m         122.0
HT tower 5 m         122.0
HT tower 8 m         122.0
HT tower 15 m        122.0
HT tower 24 m        122.0
HT tower 34 m        122.0
HT tower 49 m        122.0
HT WM                121.4
CP Centre Point      114.0
CP UK mf             113.0
CP BSE 40            114.0
CP FRG mf            114.0
                     ...  
BSE170                36.0
ANE10                107.1
ANE20                 84.5
ANE40                 41.0
ASW10                103.9
ASW20                 75.5
ASW35                 34.0
ASW50                 10.0
ASW UK 30 m tower      9.0
ASW85                 10.0
AANE10               109.8
AANE20                84.1
AANE30                59.8
AANE40                40.5
AANE50                29.0
AANE60                27.0
AASW10 mf            100.7
AASW10 t             101.0
AASW20                78.6
AASW30 mf             51.2
AASW30 t              51.4
AASW40                26.2
AASW50 mf             12.6
AASW50 t              12.0
AASW60                 9.0
AASW70                 9.1
AASW80                11.0
AASW90                10.0
small_domain           NaN
large_domain           NaN
Name: TerrainHeight, Length: 76, dtype: float64

DOWNLOAD MAPS

In [49]:
from core.wmts import WMTSDownloader
print(WMTSDownloader.list_available_servers())
['OSM', 'GEOIMAGE', 'geolandbasemap', 'Basemap AT Terrain', 'Basemap AT Overlay', 'Basemap AT HiDPI', 'Basemap AT orthofoto30cm', 'Thunderforest Outdoors', 'KOMPASS Touristik', 'KOMPASS Summer', 'KOMPASS Winter', 'Google Satellite', 'Google Road', 'Google Terrain', 'ARCGIS World_Topo_Map', 'ARCGIS World_Street_Map', 'ARCGIS World_Shaded_Relief', 'ARCGIS World_Physical_Map', 'ARCGIS World_Imagery', 'ARCGIS NatGeo_World_Map', 'APPLE', 'MapBox', 'HERE Satellite Day', 'HERE Hybrid Day', 'OpenTopoMap', 'Bing Road', 'Bing Aerial', 'Bergfex_OEK50', 'SkiTouren Winter', 'SkiTouren Summer', 'Corine 2006', 'AuWiPot Pbez', 'AuWiPot Ebez', 'VortexMap 80m', 'SWISS Topo', 'Slope overlay', 'IGN', 'ARCGIS Ocean Basemap']
In [36]:
from core.wmts.amap_downloader import AMAPDownloader
from core.wmts import WMTSDownloader
from core.docker import from_container_path_to_win

for item_cfg in ANALYSIS_CFG.maps.get("items", []):
    path = make_path(ANALYSIS_CFG.maps.get("path", os.getcwd()))
    
    pprint(item_cfg)
    domain_name = item_cfg["domain_name"]    
    domain_shape = domains_dict[domain_name]["shape"]
    mpp = item_cfg["mpp"]
    server_name = item_cfg["server_name"]
    dl = None
    if server_name == "AMAP":
        dl = AMAPDownloader(
            resolution=mpp,
            output_path=str(path),
            jpeg_quality=80,
            n_threads=32,
            optimize_download=False,
            out_epsg=ANALYSIS_CFG.projection.epsg,
            download_shape=domain_shape,
            download_shape_epsg_code=ANALYSIS_CFG.projection.epsg,
            big_tif=True,
            parallel=True,
            prefix=None,
            tag=item_cfg["filename"],
            jpeg_extension="jpeg"
        )
    else:
        print(mpp)
        dl = WMTSDownloader(
            server_name=server_name,
            resolution=mpp,
            output_path=str(path),
            force_resolution=False,
            jpeg_quality=80,
            n_threads=32,
            optimize_download=False,
            max_attempts=10,
            out_epsg=ANALYSIS_CFG.projection.epsg,
            download_shape=domain_shape,
            download_shape_epsg_code=ANALYSIS_CFG.projection.epsg,
            big_tif=True,
            parallel=True,
            jpeg_extension="jpeg"
        )
        dl.tag_name = item_cfg["filename"]
    #if os.path.isfile(dl.jpeg_fname):
    #    continue   
        
    try:  
        print(dl.download_df.index.size, getattr(dl, "selected_zoom", None), getattr(dl, "selected_resolution", None))
        dl.run_all(create_kmz=False, create_tif=False, create_jpeg=True, create_png=False, clean_files=True, create_pyramids=False, verbose=False)
    except:
        raise
        continue
    fig, ax = plt.subplots(1, 1, figsize=(20, 20))
    nrg = np.dstack((rio.open(dl.jpeg_fname).read(1), rio.open(dl.jpeg_fname).read(2), rio.open(dl.jpeg_fname).read(3)))
    ax.imshow(nrg)
    ax.set_aspect("equal")
    IPython.display.display(fig)
    plt.close(fig)
{'domain_name': 'large_domain',
 'filename': 'google_satellite_large',
 'mpp': 15,
 'server_name': 'Google Satellite'}
15
196 13 10.357841018701048
{'domain_name': 'large_domain',
 'filename': 'shaded_relief_large',
 'mpp': 15,
 'server_name': 'ARCGIS World_Shaded_Relief'}
15
196 13 10.357841018701048
{'domain_name': 'small_domain',
 'filename': 'google_satellite_small',
 'mpp': 8,
 'server_name': 'Google Satellite'}
8
36 13 10.357814179319142
In [37]:
from wind.gramm_gral.utils import change_world_files
path = make_path(ANALYSIS_CFG.maps.get("path", os.getcwd()))    
list(path.glob("*.jpegw"))
Out[37]:
[PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/karten/google_satellite_large.jpegw'),
 PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/karten/google_satellite_small.jpegw'),
 PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/windpro/basis/karten/shaded_relief_large.jpegw')]
In [38]:
def change_world_files(maps_path):
    maps_path = make_path(maps_path)
    assert maps_path.exists()
    n = 0
    for (ext, wld) in [("tif", "tfw"), ("jpeg", "jpegw"), ("png", "pngw")]:
        for fname in maps_path.glob("*." + wld):
            picture_file = fname.parent.joinpath(fname.stem + "." + ext)
            fname_new = fname.parent.joinpath(fname.stem + "." + "jpgw")
            if not picture_file.exists():
                continue
            n += 1
            picture_file = from_container_path_to_win(str(picture_file))  
            lines = fname.open("r").readlines()[:6]
            lines.append(picture_file)
            lines = "".join(lines)            
            with fname_new.open("w") as fout:
                fout.write(lines)
    return n

change_world_files(path)
Out[38]:
3

CREATE GRAMM SIMULATION

In [61]:
SIMULATION_PATH = pathlib.Path(norm_pathname(r"T:\EWS_Projects\Gral_Gramm\askervein\gramm"))
assert SIMULATION_PATH.exists()
pprint(list(SIMULATION_PATH.glob("*")))
from wind.gramm_gral.gramm.simulation import GrammSimulation
simulation = GrammSimulation(SIMULATION_PATH, epsg_code=ANALYSIS_CFG.projection.epsg)
print(simulation.projection_name)
simulation
[PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/gramm/Computation'),
 PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/gramm/Emissions'),
 PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/gramm/Maps'),
 PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/gramm/Metfiles'),
 PosixPath('/home/ewsuser/ews_drives/t_EWS_Projects/Gral_Gramm/askervein/gramm/Settings')]
OSGB 1936 / British National Grid
Out[61]:
<GrammSimulation name='gramm' width=12.9km height=12.6km resolution=50.0m ncells=2860704 />
In [62]:
gdf.query("Type != 'Domain'")
Out[62]:
AGL geometry Type TerrainHeight
Name
RS tower 1 m 1.0 POINT (74300 820980) NaN 0.0
RS tower 3 m 3.0 POINT (74300 820980) NaN 0.0
RS tower 5 m 5.0 POINT (74300 820980) NaN 0.0
RS tower 8 m 8.0 POINT (74300 820980) NaN 0.0
RS tower 15 m 15.0 POINT (74300 820980) NaN 0.0
RS tower 24 m 24.0 POINT (74300 820980) NaN 0.0
RS tower 34 m 34.0 POINT (74300 820980) NaN 0.0
RS tower 49 m 49.0 POINT (74300 820980) NaN 0.0
Cemetary BM -10000.0 POINT (73380 821940) NaN 14.2
Milestone -10000.0 POINT (75284 821714) NaN 7.0
TV tower -10000.0 POINT (73590 821711) NaN 6.3
Base station -10000.0 POINT (74846 823306) NaN 9.0
Pacman lake -10000.0 POINT (76207 823587) NaN 38.0
Askernish House -10000.0 POINT (73654 823845) NaN 4.0
HT Hill Top 10.0 POINT (75383 823737) NaN 122.0
HT 10 m mf 10.0 POINT (75387 823735) NaN 122.0
HT 10 m t 10.0 POINT (75381 823745) NaN 122.0
HT tower 1 m 1.0 POINT (75381 823753) NaN 122.0
HT tower 3 m 3.0 POINT (75381 823753) NaN 122.0
HT tower 5 m 5.0 POINT (75381 823753) NaN 122.0
HT tower 8 m 8.0 POINT (75381 823753) NaN 122.0
HT tower 15 m 15.0 POINT (75381 823753) NaN 122.0
HT tower 24 m 24.0 POINT (75381 823753) NaN 122.0
HT tower 34 m 34.0 POINT (75381 823753) NaN 122.0
HT tower 49 m 49.0 POINT (75381 823753) NaN 122.0
HT WM 10.0 POINT (75360 823763) NaN 121.4
CP Centre Point 10.0 POINT (75678 823465) NaN 114.0
CP UK mf 10.0 POINT (75686 823457) NaN 113.0
CP BSE 40 10.0 POINT (75680 823464) NaN 114.0
CP FRG mf 10.0 POINT (75676 823466) NaN 114.0
... ... ... ... ...
BSE130 10.0 POINT (76338 822860) NaN 31.0
BSE150 10.0 POINT (76490 822723) NaN 31.0
BSE170 10.0 POINT (76636 822585) NaN 36.0
ANE10 10.0 POINT (75454 823812) NaN 107.1
ANE20 10.0 POINT (75523 823884) NaN 84.5
ANE40 10.0 POINT (75661 824017) NaN 41.0
ASW10 10.0 POINT (75319 823667) NaN 103.9
ASW20 10.0 POINT (75257 823604) NaN 75.5
ASW35 10.0 POINT (75162 823498) NaN 34.0
ASW50 10.0 POINT (75050 823378) NaN 10.0
ASW UK 30 m tower 30.0 POINT (74970 823295) NaN 9.0
ASW85 10.0 POINT (74813 823122) NaN 10.0
AANE10 10.0 POINT (75746 823540) NaN 109.8
AANE20 10.0 POINT (75807 823610) NaN 84.1
AANE30 10.0 POINT (75871 823675) NaN 59.8
AANE40 10.0 POINT (75938 823745) NaN 40.5
AANE50 10.0 POINT (76006 823813) NaN 29.0
AANE60 10.0 POINT (76073 823886) NaN 27.0
AASW10 mf 10.0 POINT (75613 823396) NaN 100.7
AASW10 t 10.0 POINT (75623 823388) NaN 101.0
AASW20 10.0 POINT (75553 823320) NaN 78.6
AASW30 mf 10.0 POINT (75483 823253) NaN 51.2
AASW30 t 10.0 POINT (75493 823244) NaN 51.4
AASW40 10.0 POINT (75417 823174) NaN 26.2
AASW50 mf 10.0 POINT (75341 823107) NaN 12.6
AASW50 t 10.0 POINT (75352 823100) NaN 12.0
AASW60 10.0 POINT (75274 823038) NaN 9.0
AASW70 10.0 POINT (75208 822968) NaN 9.1
AASW80 10.0 POINT (75140 822895) NaN 11.0
AASW90 10.0 POINT (75069 822820) NaN 10.0

74 rows × 4 columns

In [79]:
print(simulation.write_receptors(gdf.query("Type != 'Domain'").assign(Height=10, Diameter=10), height_name="Height", diameter_name="Diameter"))
74
1,74300,820980,10,RS_tower_1_m,10
2,74300,820980,10,RS_tower_3_m,10
3,74300,820980,10,RS_tower_5_m,10
4,74300,820980,10,RS_tower_8_m,10
5,74300,820980,10,RS_tower_15_m,10
6,74300,820980,10,RS_tower_24_m,10
7,74300,820980,10,RS_tower_34_m,10
8,74300,820980,10,RS_tower_49_m,10
9,73380,821940,10,Cemetary_BM,10
10,75284,821714,10,Milestone,10
11,73590,821711,10,TV_tower,10
12,74846,823306,10,Base_station,10
13,76207,823587,10,Pacman_lake,10
14,73654,823845,10,Askernish_House,10
15,75383,823737,10,HT_Hill_Top,10
16,75387,823735,10,HT_10_m_mf,10
17,75381,823745,10,HT_10_m_t,10
18,75381,823753,10,HT_tower_1_m,10
19,75381,823753,10,HT_tower_3_m,10
20,75381,823753,10,HT_tower_5_m,10
21,75381,823753,10,HT_tower_8_m,10
22,75381,823753,10,HT_tower_15_m,10
23,75381,823753,10,HT_tower_24_m,10
24,75381,823753,10,HT_tower_34_m,10
25,75381,823753,10,HT_tower_49_m,10
26,75360,823763,10,HT_WM,10
27,75678,823465,10,CP_Centre_Point,10
28,75686,823457,10,CP_UK_mf,10
29,75680,823464,10,CP_BSE_40,10
30,75676,823466,10,CP_FRG_mf,10
31,75688,823493,10,CP_FRG_17_m_tower,10
32,75313,823810,10,BNW10,10
33,75243,823875,10,BNW20,10
34,75458,823671,10,BSE10,10
35,75528,823603,10,BSE20,10
36,75606,823535,10,BSE30,10
37,75680,823465,10,BSE40,10
38,75754,823397,10,BSE50,10
39,75833,823324,10,BSE60,10
40,75905,823260,10,BSE70,10
41,75982,823193,10,BSE80,10
42,76046,823130,10,BSE90,10
43,76120,823063,10,BSE100,10
44,76195,822997,10,BSE110,10
45,76338,822860,10,BSE130,10
46,76490,822723,10,BSE150,10
47,76636,822585,10,BSE170,10
48,75454,823812,10,ANE10,10
49,75523,823884,10,ANE20,10
50,75661,824017,10,ANE40,10
51,75319,823667,10,ASW10,10
52,75257,823604,10,ASW20,10
53,75162,823498,10,ASW35,10
54,75050,823378,10,ASW50,10
55,74970,823295,10,ASW_UK_30_m_tower,10
56,74813,823122,10,ASW85,10
57,75746,823540,10,AANE10,10
58,75807,823610,10,AANE20,10
59,75871,823675,10,AANE30,10
60,75938,823745,10,AANE40,10
61,76006,823813,10,AANE50,10
62,76073,823886,10,AANE60,10
63,75613,823396,10,AASW10_mf,10
64,75623,823388,10,AASW10_t,10
65,75553,823320,10,AASW20,10
66,75483,823253,10,AASW30_mf,10
67,75493,823244,10,AASW30_t,10
68,75417,823174,10,AASW40,10
69,75341,823107,10,AASW50_mf,10
70,75352,823100,10,AASW50_t,10
71,75274,823038,10,AASW60,10
72,75208,822968,10,AASW70,10
73,75140,822895,10,AASW80,10
74,75069,822820,10,AASW90,10

In [64]:
simulation.mesh.read_mesh()
simulation.mesh.write_terrain_data(simulation.project_path.joinpath("terrain_data.vti"))
simulation.mesh.write_vtk_grid(simulation.project_path.joinpath("mesh.vts"))
print(simulation.number_of_cells, simulation.domain_width, simulation.domain_height)
print(simulation.receptors_gdf.distance(simulation.domain.values[0].exterior).min())
print(simulation.latitude)
print(simulation.update_config_files_with_latitude())
/usr/local/lib/python3.6/site-packages/vtk/util/numpy_support.py:137: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.
  assert not numpy.issubdtype(z.dtype, complex), \
2860704 12900.0 12600.0
4730.0
57.17820164744275
None
In [65]:
simulation.read_computation_info()
print(simulation.uses_meteopgt)
print(simulation.z0)
print(simulation.n_smoothed_cells)
print(simulation.write_steady_state_file)
True
0.2
20
False
In [69]:
simulation.n_smoothed_cells
In [70]:
rasters = simulation.topography_raster, simulation.land_use_raster, simulation.tpi_slope_min_raster, simulation.landform_class_raster
n = int(len([_ for _ in rasters if _ is not None]) / 2)
fig, axL = plt.subplots(n, 2, figsize=(8 * 2, n * (8*simulation.aspect_ratio + 0.5) ))
axL = np.array(axL).flatten()
simulation._rasters_dict = dict()
bounds = simulation.bounds
extent = simulation.extent
for ax, raster in zip(axL, rasters):
    cmap = "gist_earth" if raster.name == "topography" else "viridis"
    norm = None
    if raster.name == "landform_class":
        import matplotlib as mpl
        from core.eplot.pandas_integration import EchartsAxis
        cmap = mpl.colors.ListedColormap(EchartsAxis.create_palette("tab10_10"))
        norm = mpl.colors.BoundaryNorm(np.arange(1, 11, 1), cmap.N) 
    im = ax.imshow(raster.data, cmap=cmap, norm=norm, extent=raster.extent)
    
    plt.colorbar(im, ax=ax, orientation="horizontal").set_label(raster.name)
    if True:
        simulation.n_smoothed_cells = 20
        ax.plot(
            simulation.receptors_gdf.X, 
            simulation.receptors_gdf.Y, 
            mfc="w", color="k", marker="o", ls="none"
        )
        smoothing_distance = simulation.dx * simulation.n_smoothed_cells
        style = dict(color="white", ls="--")
        ax.axvline(bounds[0] + smoothing_distance, **style)
        ax.axvline(bounds[2] - smoothing_distance, **style)
        ax.axhline(bounds[1] + smoothing_distance, **style)
        ax.axhline(bounds[3] - smoothing_distance, **style)
    if True:
        distance = 10000
        style = dict(color="red", ls="--")
        ax.axvline(bounds[0] + distance, **style)
        ax.axvline(bounds[2] - distance, **style)
        ax.axhline(bounds[1] + distance, **style)
        ax.axhline(bounds[3] - distance, **style)
        
    ax.set(xlim=extent[:2], ylim=extent[2:], aspect="equal", title=raster.name)
fig.tight_layout();
In [71]:
simulation.write_meteopgt(n_sectors=36, wind_speeds=[5], stability_classes=[4], height=10)
simulation.test_cases_df
Out[71]:
WindDirectionSector WindSpeedClass StabilityClass Frequency WindDirection
testcase
00001 0.0 5.0 4.0 0.0036 0.0
00002 1.0 5.0 4.0 0.0036 10.0
00003 2.0 5.0 4.0 0.0036 20.0
00004 3.0 5.0 4.0 0.0036 30.0
00005 4.0 5.0 4.0 0.0036 40.0
00006 5.0 5.0 4.0 0.0036 50.0
00007 6.0 5.0 4.0 0.0036 60.0
00008 7.0 5.0 4.0 0.0036 70.0
00009 8.0 5.0 4.0 0.0036 80.0
00010 9.0 5.0 4.0 0.0036 90.0
00011 10.0 5.0 4.0 0.0036 100.0
00012 11.0 5.0 4.0 0.0036 110.0
00013 12.0 5.0 4.0 0.0036 120.0
00014 13.0 5.0 4.0 0.0036 130.0
00015 14.0 5.0 4.0 0.0036 140.0
00016 15.0 5.0 4.0 0.0036 150.0
00017 16.0 5.0 4.0 0.0036 160.0
00018 17.0 5.0 4.0 0.0036 170.0
00019 18.0 5.0 4.0 0.0036 180.0
00020 19.0 5.0 4.0 0.0036 190.0
00021 20.0 5.0 4.0 0.0036 200.0
00022 21.0 5.0 4.0 0.0036 210.0
00023 22.0 5.0 4.0 0.0036 220.0
00024 23.0 5.0 4.0 0.0036 230.0
00025 24.0 5.0 4.0 0.0036 240.0
00026 25.0 5.0 4.0 0.0036 250.0
00027 26.0 5.0 4.0 0.0036 260.0
00028 27.0 5.0 4.0 0.0036 270.0
00029 28.0 5.0 4.0 0.0036 280.0
00030 29.0 5.0 4.0 0.0036 290.0
00031 30.0 5.0 4.0 0.0036 300.0
00032 31.0 5.0 4.0 0.0036 310.0
00033 32.0 5.0 4.0 0.0036 320.0
00034 33.0 5.0 4.0 0.0036 330.0
00035 34.0 5.0 4.0 0.0036 340.0
00036 35.0 5.0 4.0 0.0036 350.0
In [77]:
simulation.number_of_cells
Out[77]:
2860704
In [72]:
simulation.clean_computation_files()
Out[72]:
15
In [76]:
simulation.mesh.read_mesh()
simulation.clean_computation_files()
simulation.clean_jobs_files()
simulation.update_config_files_with_latitude()
simulation.set_use_flat(False)

for i, l in enumerate(simulation.computation_path.joinpath("IIN.dat").open("r").readlines()):
    print(f"{i}:{l.strip()}")

test_cases_df = simulation.test_cases_df
print(test_cases_df.index.size)
jobs_df = simulation.write_job_files(n_procs=12, n_threads=1)
#jobs_df = simulation.write_job_files(n_procs=5, n_threads=5, query=None)
print(test_cases_df.index.size, jobs_df.index.size, jobs_df.index[0])
jobs_df
0:COMPUTATION DATE         (YYMMDD)             :  080210               !No influence when using meteopgt.all
1:BEGINNING OF COMPUTATION   (hhmm)             :  0000                 !No influence when using meteopgt.all
2:MAXIMUM ALLOWED TIME STEP DT  [ s ]           :  10                   !Range 1-100
3:MODELLING TIME (FOR VALUES >1(s) AND <1[%])   :  3600                 !Range 0.01-1% or 2-infinite sec.
4:BUFFERING AFTER TIMESTEPS                     :  3600                 !No influence when using meteopgt.all
5:MAX. PERMISSIBLE W-DEVIATION ABOVE < 1 [mm/s] :  0.01                 !Not used
6:RELATIVE INITIAL HUMIDITY  [ % ] GT.0         :  20                   !No influence when using meteopgt.all
7:HEIGHT OF THE LOWEST COMPUTATION HEIGHT [ m ] :  330.                 !Not used
8:AIR TEMPERATURE AT GROUND  [ K ]              :  280.0                !No influence when using meteopgt.all
9:GRADIENT OF THE POTENTIAL TEMPERATURE [K/100m]:  -0.0065              !Not used
10:NEUTRAL LAYERING UP TO THE HEIGHT ABOVE GROUND:  4000                 !Not used
11:SURFACE TEMPERATURE [ K ]                     :  280.0                !No influence when using meteopgt.all
12:TEMPERATURE OF THE SOIL IN 1 M DEPTH  [ K ]   :  280.0                !No influence when using meteopgt.all
13:LATITUDE                                      :  57                   !Range -90 - +90 deg.
14:UPDATE OF RADIATION (TIMESTEPS)               :  300                  !No influence when using meteopgt.all
15:DEBUG LEVEL 0 none, 3 highest                 :  0                    !not used
16:COMPUTE  U V W PN T GW FO            YES=1    :  1111000              !GW and FO not used
17:COMPUTE  BR PR QU PSI TE TB STR      YES=1    :  1111111              !BR, PSI not used, TE not recommended
18:DIAGNOSTIC INITIAL CONDITIONS        YES=1    :  1                    !Not used
19:EXPLICIT/IMPLICIT TIME INTEGRATION  I=1/E=0   :  1                    !Not used
20:RELAXATION VELOCITY                           :  0.1                  !Range 0.01-1.0
21:RELAXATION SCALARS                            :  0.1                  !Range 0.01-1.0
22:Force catab flows -40/-25/-15 W/m² AKLA 7/6/5 :  0                    !1=ON suitable with steady state meteo.all
23:NESTING (=1) OR LARGE SCALE FORCING    (=0)   :  0                    !Not used
24:BOUNDARY CONDITION (1,2,3,4,5,6)              :  6                    !1=homogenous Von Neumann, 6=Nudging (standard)
25:NON-STEADYSTATE/FORCING W. 'METEO.ALL' YES=1  :  0                    !Should not be changed
26:NON-STEADYSTATE/FORCING W. IFU-MM5 DATA YES=1 :  0                    !Not used
27:GRAMM IN GRAMM NESTING YES=1                  :  0                    !Not used
28:Flowfield Output Format                       :  0                    !0=binary stream +header, 3=binary stream
36
36 36 00001
Out[76]:
WindDirectionSector WindSpeedClass StabilityClass Frequency WindDirection processor
testcase
00001 0.0 5.0 4.0 0.0036 0.0 0
00013 12.0 5.0 4.0 0.0036 120.0 0
00025 24.0 5.0 4.0 0.0036 240.0 0
00002 1.0 5.0 4.0 0.0036 10.0 1
00014 13.0 5.0 4.0 0.0036 130.0 1
00026 25.0 5.0 4.0 0.0036 250.0 1
00003 2.0 5.0 4.0 0.0036 20.0 2
00015 14.0 5.0 4.0 0.0036 140.0 2
00027 26.0 5.0 4.0 0.0036 260.0 2
00004 3.0 5.0 4.0 0.0036 30.0 3
00016 15.0 5.0 4.0 0.0036 150.0 3
00028 27.0 5.0 4.0 0.0036 270.0 3
00005 4.0 5.0 4.0 0.0036 40.0 4
00017 16.0 5.0 4.0 0.0036 160.0 4
00029 28.0 5.0 4.0 0.0036 280.0 4
00006 5.0 5.0 4.0 0.0036 50.0 5
00018 17.0 5.0 4.0 0.0036 170.0 5
00030 29.0 5.0 4.0 0.0036 290.0 5
00007 6.0 5.0 4.0 0.0036 60.0 6
00019 18.0 5.0 4.0 0.0036 180.0 6
00031 30.0 5.0 4.0 0.0036 300.0 6
00008 7.0 5.0 4.0 0.0036 70.0 7
00020 19.0 5.0 4.0 0.0036 190.0 7
00032 31.0 5.0 4.0 0.0036 310.0 7
00009 8.0 5.0 4.0 0.0036 80.0 8
00021 20.0 5.0 4.0 0.0036 200.0 8
00033 32.0 5.0 4.0 0.0036 320.0 8
00010 9.0 5.0 4.0 0.0036 90.0 9
00022 21.0 5.0 4.0 0.0036 210.0 9
00034 33.0 5.0 4.0 0.0036 330.0 9
00011 10.0 5.0 4.0 0.0036 100.0 10
00023 22.0 5.0 4.0 0.0036 220.0 10
00035 34.0 5.0 4.0 0.0036 340.0 10
00012 11.0 5.0 4.0 0.0036 110.0 11
00024 23.0 5.0 4.0 0.0036 230.0 11
00036 35.0 5.0 4.0 0.0036 350.0 11
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [82]:
display_end_of_script(raise_error=False, tz='Europe/Vienna')
END OF THE SCRIPT
2019-11-08 10:17:34.865124+01:00
In [83]:
IPython.display.HTML('''<script>
code_show = true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
</script>
<button id="__toggler__" style="position: fixed;right: 5px;top: 5px;z-index: 1000;" type="button" onclick="code_toggle()"><b>TOGGLE CODE CELLS</b></button>''')
Out[83]:
In [ ]:
nbconvert_html('00_setup_windpro_project.ipynb', toc=True)
In [1]:
!jupyter nbconvert --to=html_toc --ExtractOutputPreprocessor.extract_output_types=[] --ExecutePreprocessor.store_widget_state=True '00_setup_windpro_project.ipynb'
[NbConvertApp] Converting notebook 00_setup_windpro_project.ipynb to html_toc
[NbConvertApp] Writing 12515195 bytes to 00_setup_windpro_project.html
In [ ]: